OpenFlights.org (https://www.openflights.org/data.html) が提供している航空網のデータセットを対象としたネットワーク分析を行う。 航空網ネットワークの性質を調べたり、航空網のつながり方を用いた空港の分類や重要な空港の発見を行う。
import pandas as pd
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from networkx.algorithms.community import greedy_modularity_communities
from collections import Counter
import random
import pickle
%matplotlib inline
airports-extended.dat, routes.datの中身を確認する。¶中身のテキストを確認できたらOK!
airports-extended.dat, routes.dat をpandasのdataframeとして読み込む。¶さらに、https://www.openflights.org/data.html を参考に、各列に適切な名前をつける。
(ヒント)
airports-extended.dat はカンマ(,)で各列が区切られた表形式のテキストデータになっている。CSV(Comma-Separated Values) という。read_csvメソッドを用いる。airports = pd.read_csv("airports-extended.dat", # ここはファイルの置いてある場所に合わせて変える
names=[
"Airport ID",
"Name",
"City",
"Country",
"IATA",
"ICAO",
"Latitude",
"Longitude",
"Altitude",
"Timezone",
"DST",
"Tz database time zone",
"Type",
"Source"
])
airports.head()
| Airport ID | Name | City | Country | IATA | ICAO | Latitude | Longitude | Altitude | Timezone | DST | Tz database time zone | Type | Source | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | Goroka Airport | Goroka | Papua New Guinea | GKA | AYGA | -6.081690 | 145.391998 | 5282 | 10 | U | Pacific/Port_Moresby | airport | OurAirports |
| 1 | 2 | Madang Airport | Madang | Papua New Guinea | MAG | AYMD | -5.207080 | 145.789001 | 20 | 10 | U | Pacific/Port_Moresby | airport | OurAirports |
| 2 | 3 | Mount Hagen Kagamuga Airport | Mount Hagen | Papua New Guinea | HGU | AYMH | -5.826790 | 144.296005 | 5388 | 10 | U | Pacific/Port_Moresby | airport | OurAirports |
| 3 | 4 | Nadzab Airport | Nadzab | Papua New Guinea | LAE | AYNZ | -6.569803 | 146.725977 | 239 | 10 | U | Pacific/Port_Moresby | airport | OurAirports |
| 4 | 5 | Port Moresby Jacksons International Airport | Port Moresby | Papua New Guinea | POM | AYPY | -9.443380 | 147.220001 | 146 | 10 | U | Pacific/Port_Moresby | airport | OurAirports |
routes = pd.read_csv("routes.dat", # ここはファイルの置いてある場所に合わせて変える
names=[
"Airline",
"Airline ID",
"Source airport",
"Source airport ID",
"Destination airport",
"Destination airport ID",
"Codeshare",
"Stops",
"Equipment"
])
routes.head()
| Airline | Airline ID | Source airport | Source airport ID | Destination airport | Destination airport ID | Codeshare | Stops | Equipment | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 2B | 410 | AER | 2965 | KZN | 2990 | NaN | 0 | CR2 |
| 1 | 2B | 410 | ASF | 2966 | KZN | 2990 | NaN | 0 | CR2 |
| 2 | 2B | 410 | ASF | 2966 | MRV | 2962 | NaN | 0 | CR2 |
| 3 | 2B | 410 | CEK | 2968 | KZN | 2990 | NaN | 0 | CR2 |
| 4 | 2B | 410 | CEK | 2968 | OVB | 4078 | NaN | 0 | CR2 |
(ヒント)
networkxライブラリのGraph形式のオブジェクトを作成する」こと。pandasライブラリのDataFramenetworkxライブラリのGraphfrom_pandas_edgelist
# ネットワークの作成
G = nx.from_pandas_edgelist(df=routes, source="Source airport", target="Destination airport", create_using=nx.Graph)
(ヒント)
draw を用いる。
pos 情報を dictionary(辞書)形式 で与える必要があるとわかる。airports-extended.datから読み込んだデータ)から緯度・経度情報を抽出し、nx.draw の pos 情報に与えてあげる。# 緯度・経度の情報を airports データフレームから取得する
pos = {data["IATA"]: (data["Longitude"], data["Latitude"]) for index, data in airports.iterrows()}
# posの中身を5個だけ見てみる
for k, v in list(pos.items())[:5]:
print(k, v)
GKA (145.391998291, -6.081689834590001) MAG (145.789001465, -5.20707988739) HGU (144.29600524902344, -5.826789855957031) LAE (146.725977, -6.569803) POM (147.22000122070312, -9.44338035583496)
# 緯度・経度情報を用いて描画
fig = plt.figure()
ax = fig.add_subplot(111)
nx.draw(G, pos, ax)
plt.show()
--------------------------------------------------------------------------- KeyError Traceback (most recent call last) ~/.anyenv/envs/pyenv/versions/miniforge3-4.9.2/lib/python3.9/site-packages/networkx/drawing/nx_pylab.py in draw_networkx_nodes(G, pos, nodelist, node_size, node_color, node_shape, alpha, cmap, vmin, vmax, ax, linewidths, edgecolors, label) 455 try: --> 456 xy = np.asarray([pos[v] for v in nodelist]) 457 except KeyError as e: ~/.anyenv/envs/pyenv/versions/miniforge3-4.9.2/lib/python3.9/site-packages/networkx/drawing/nx_pylab.py in <listcomp>(.0) 455 try: --> 456 xy = np.asarray([pos[v] for v in nodelist]) 457 except KeyError as e: KeyError: 'INC' The above exception was the direct cause of the following exception: NetworkXError Traceback (most recent call last) <ipython-input-9-2cfa2fc01630> in <module> 3 ax = fig.add_subplot(111) 4 ----> 5 nx.draw(G, pos, ax) 6 plt.show() ~/.anyenv/envs/pyenv/versions/miniforge3-4.9.2/lib/python3.9/site-packages/networkx/drawing/nx_pylab.py in draw(G, pos, ax, **kwds) 121 kwds["with_labels"] = "labels" in kwds 122 --> 123 draw_networkx(G, pos=pos, ax=ax, **kwds) 124 ax.set_axis_off() 125 plt.draw_if_interactive() ~/.anyenv/envs/pyenv/versions/miniforge3-4.9.2/lib/python3.9/site-packages/networkx/drawing/nx_pylab.py in draw_networkx(G, pos, arrows, with_labels, **kwds) 333 pos = nx.drawing.spring_layout(G) # default to spring layout 334 --> 335 draw_networkx_nodes(G, pos, **node_kwds) 336 draw_networkx_edges(G, pos, arrows=arrows, **edge_kwds) 337 if with_labels: ~/.anyenv/envs/pyenv/versions/miniforge3-4.9.2/lib/python3.9/site-packages/networkx/drawing/nx_pylab.py in draw_networkx_nodes(G, pos, nodelist, node_size, node_color, node_shape, alpha, cmap, vmin, vmax, ax, linewidths, edgecolors, label) 456 xy = np.asarray([pos[v] for v in nodelist]) 457 except KeyError as e: --> 458 raise nx.NetworkXError(f"Node {e} has no position.") from e 459 except ValueError as e: 460 raise nx.NetworkXError("Bad value in node positions.") from e NetworkXError: Node 'INC' has no position.
エラーがでる。おそらく原因は、各空港データ (airports-extended.dat)にすべての空港の緯度・経度が含まれていないせい。
仕方がないので、今回の描画はそこに含まれる空港のみを抽出して行う。
# ネットワーク G から、緯度・経度情報がある空港ノードのみを抽出したサブグラフを作成
G_plot = G.subgraph(pos.keys())
# あらためて描画
fig = plt.figure()
ax = fig.add_subplot(111)
nx.draw(G_plot, pos, ax)
plt.show()
# 見づらいので、いくつかパラメータを与えて少しだけ整える
fig = plt.figure(figsize=(15,7))
ax = fig.add_subplot(111)
nx.draw(G_plot, pos, ax, node_size=30, style=":")
plt.show()
ちなみに、ネットワーク構造の可視化は Gephi や Cytoscape といったソフトウェアを用いて行うのもよい
# 基本統計量
print(nx.info(G))
Name: Type: Graph Number of nodes: 3425 Number of edges: 19257 Average degree: 11.2450
縦軸・横軸は対数軸でプロットする。
(ヒント)
# 次数分布
## 次数の集計
degree_list = [degree for name, degree in G.degree(G)]
degree_count = Counter(degree_list)
## 描画
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)
ax.scatter(degree_count.keys(), degree_count.values())
ax.set_xlabel('degree')
ax.set_ylabel('frequency')
ax.set_xscale('log')
ax.set_yscale('log')
あるネットワークのうち、任意のノード間にパス(経路)が存在するサブグラフを接続成分(connected components)という。特に、ネットワーク内で最大の接続成分を最大接続成分(largest connected component)という。
(例)以下のネットワークには、接続成分が3つある。
# 各接続成分の抽出
nx.connected_components(G)
<generator object connected_components at 0x13e711430>
connected_componentsのドキュメントを読むと、Returnsの部分に
A generator of sets of nodes, one for each component of G.
と書いてある。
ここから、connected_componentsの返り値は「各接続成分のノードセットのgenerator」であることがわかる。
generator型オブジェクトはfor文で中身を取り出すことができるほか、list()を適用することでlist型に変換することができる。
# 各接続成分のノード数をチェック
for cc in nx.connected_components(G):
print(len(cc))
3397 2 4 2 4 2 4 10
一つ目のconnected componentsが最大であるため、それを抽出する。
cc_list = list(nx.connected_components(G)) # 各connected componentsのノードセットをlist()で配列化する。
lcc_node_list = cc_list[0] # 最大接続成分(largest connected components)のノードセットを取り出す。
G_lcc = G.subgraph(lcc_node_list) # ノードセットを指定して、最大接続成分のサブグラフを元のネットワークから抽出する。
ネットワークから最大接続成分を抽出することで、平均パス長(ネットワークに含まれるノードの全ペアの最短経路長を平均した値)や、直径(ネットワーク内の最大の最短経路長)といった値を求められるようになる。これらは、ネットワークの性質を示す重要な値である。
今回の航空網ネットワークではほとんどのノードが最大接続成分に含まれるため、最大接続成分の性質を調べることでもネットワーク全体の性質を知ることができた、といえるだろう。
# 直径 (けっこう実行に時間がかかる)
nx.diameter(G_lcc)
13
# 平均パス長 (けっこう実行に時間がかかる)
nx.average_shortest_path_length(G_lcc)
4.103241167898093
あるネットワーク固有の性質を調べるためのアプローチとして、null model(帰無仮説ネットワーク)と比較する方法がある。
「現実のネットワーク」と、「null model (= あるランダムなネットワーク生成アルゴリズムによって、ノード数や次数などの基本的な性質が同じになるように再現したネットワーク)」を比較する。このとき、ネットワーク特徴の指標(直径、平均パス長、ネットワーク・モチーフ、etc...)に有意な差がみられる場合、その部分に現実のネットワークの構造を規定する制約が反映されているのではないか?と考えられる。ひいては、現実世界のシステムを理解する手がかりになる。
(ヒント)
n は最大接続成分のノード数p は最大接続成分の密度seed は与えても与えなくてもいいが、適当な値を与えておくようにすると再現ができるためよい。n = len(G_lcc.nodes)
p = nx.density(G_lcc)
seed = 2021
G_rand = nx.erdos_renyi_graph(n, p, seed)
print("元データのノード数:", len(G_lcc.nodes))
print("ランダムのノード数:", len(G_rand.nodes))
元データのノード数: 3397 ランダムのノード数: 3397
print("元データのエッジ数", len(G_lcc.edges))
print("ランダムのエッジ数", len(G_rand.edges))
元データのエッジ数 19231 ランダムのエッジ数 19165
# ランダムグラフの直径
nx.diameter(G_rand)
6
# ランダムグラフの平均パス長
nx.average_shortest_path_length(G_rand)
3.625276303868202
ノード数・エッジ数を再現したランダムグラフよりも、実際の航空網ネットワークのほうが直径・平均パス長ともに大きい。このような差を生み出す原因には、ランダムグラフで考慮されていなかった地理的要因や政治的要因などがあると考えられる。
ランダムグラフの生成アルゴリズムにそのような要因を追加して、再び現実の航空網と比較し、指標をよりよく再現できているか?を確かめることで、現実の航空網ネットワークを生み出すメカニズムをよりよく理解できる。
(注:本当はランダムグラフを複数サンプリングしたり、ランダム化した場合の値を解析的に求めたりして、より厳密にやる必要があるかもしれない。)
# 次数中心性
degree_centralities = nx.degree_centrality(G)
# ネットワークのデータに追記する
nx.set_node_attributes(G, degree_centralities, "degree_centrality")
# 各空港データにも結合する
degree_centralities_df = pd.DataFrame(degree_centralities.items(), columns=["IATA", "degree_centrality"])
airports = pd.merge(airports, degree_centralities_df, on="IATA", how="left")
# 次数中心性の上位10個がどんな空港なのか調べる
airports.sort_values("degree_centrality", ascending=False).iloc[:10]
| Airport ID | Name | City | Country | IATA | ICAO | Latitude | Longitude | Altitude | Timezone | DST | Tz database time zone | Type | Source | degree_centrality | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 574 | 580 | Amsterdam Airport Schiphol | Amsterdam | Netherlands | AMS | EHAM | 52.308601 | 4.763890 | -11 | 1 | E | Europe/Amsterdam | airport | OurAirports | 0.072430 |
| 336 | 340 | Frankfurt am Main Airport | Frankfurt | Germany | FRA | EDDF | 50.033333 | 8.570556 | 364 | 1 | E | Europe/Berlin | airport | OurAirports | 0.071262 |
| 1350 | 1382 | Charles de Gaulle International Airport | Paris | France | CDG | LFPG | 49.012798 | 2.550000 | 392 | 1 | E | Europe/Paris | airport | OurAirports | 0.070093 |
| 12253 | 13696 | Istanbul Airport | Istanbul | Turkey | IST | LTFM | 41.275278 | 28.751944 | 325 | 3 | E | \N | airport | OurAirports | 0.068925 |
| 3494 | 3682 | Hartsfield Jackson Atlanta International Airport | Atlanta | United States | ATL | KATL | 33.636700 | -84.428101 | 1026 | -5 | A | America/New_York | airport | OurAirports | 0.063376 |
| 3181 | 3364 | Beijing Capital International Airport | Beijing | China | PEK | ZBAA | 40.080101 | 116.584999 | 116 | 8 | U | Asia/Shanghai | airport | OurAirports | 0.060456 |
| 3642 | 3830 | Chicago O'Hare International Airport | Chicago | United States | ORD | KORD | 41.978600 | -87.904800 | 672 | -6 | A | America/Chicago | airport | OurAirports | 0.060164 |
| 342 | 346 | Munich Airport | Munich | Germany | MUC | EDDM | 48.353802 | 11.786100 | 1487 | 1 | E | Europe/Berlin | airport | OurAirports | 0.056075 |
| 3830 | 4029 | Domodedovo International Airport | Moscow | Russia | DME | UUDD | 55.408798 | 37.906300 | 588 | 3 | N | Europe/Moscow | airport | OurAirports | 0.055491 |
| 2106 | 2188 | Dubai International Airport | Dubai | United Arab Emirates | DXB | OMDB | 25.252800 | 55.364399 | 62 | 4 | U | Asia/Dubai | airport | OurAirports | 0.054907 |
# 例によって緯度・経度のデータ(pos)がある空港データに絞る
G_plot = G.subgraph(pos.keys())
# centralityの値を描画するノードのサイズに利用するために抽出。 1000は適当(描画の見栄えがよくなるように)
node_size_list = [node_data["degree_centrality"] * 1000 for node_name, node_data in G_plot.nodes(data=True)]
# 描画
fig = plt.figure(figsize=(15,7))
ax = fig.add_subplot(111)
nx.draw(G_plot, pos, ax, node_size=node_size_list, style=":")
plt.show()
# 媒介中心性 少し実行に時間がかかる!
betweenness_centralities = nx.betweenness_centrality(G)
# ネットワークのデータに追記する
nx.set_node_attributes(G, betweenness_centralities, "betweenness_centrality")
# 各空港データにも結合する
betweenness_centralities_df = pd.DataFrame(betweenness_centralities.items(), columns=["IATA", "betweenness_centrality"])
airports = pd.merge(airports, betweenness_centralities_df, on="IATA", how="left")
# 媒介中心性の上位10個がどんな空港なのか調べる
airports.sort_values("betweenness_centrality", ascending=False).iloc[:10]
| Airport ID | Name | City | Country | IATA | ICAO | Latitude | Longitude | Altitude | Timezone | DST | Tz database time zone | Type | Source | degree_centrality | betweenness_centrality | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 3586 | 3774 | Ted Stevens Anchorage International Airport | Anchorage | United States | ANC | PANC | 61.174400 | -149.996002 | 152 | -9 | A | America/Anchorage | airport | OurAirports | 0.009930 | 0.072483 |
| 3296 | 3484 | Los Angeles International Airport | Los Angeles | United States | LAX | KLAX | 33.942501 | -118.407997 | 125 | -8 | A | America/Los_Angeles | airport | OurAirports | 0.043516 | 0.065208 |
| 1350 | 1382 | Charles de Gaulle International Airport | Paris | France | CDG | LFPG | 49.012798 | 2.550000 | 392 | 1 | E | Europe/Paris | airport | OurAirports | 0.070093 | 0.061936 |
| 2106 | 2188 | Dubai International Airport | Dubai | United Arab Emirates | DXB | OMDB | 25.252800 | 55.364399 | 62 | 4 | U | Asia/Dubai | airport | OurAirports | 0.054907 | 0.056488 |
| 336 | 340 | Frankfurt am Main Airport | Frankfurt | Germany | FRA | EDDF | 50.033333 | 8.570556 | 364 | 1 | E | Europe/Berlin | airport | OurAirports | 0.071262 | 0.051086 |
| 3181 | 3364 | Beijing Capital International Airport | Beijing | China | PEK | ZBAA | 40.080101 | 116.584999 | 116 | 8 | U | Asia/Shanghai | airport | OurAirports | 0.060456 | 0.049569 |
| 574 | 580 | Amsterdam Airport Schiphol | Amsterdam | Netherlands | AMS | EHAM | 52.308601 | 4.763890 | -11 | 1 | E | Europe/Amsterdam | airport | OurAirports | 0.072430 | 0.049262 |
| 3389 | 3577 | Seattle Tacoma International Airport | Seattle | United States | SEA | KSEA | 47.449001 | -122.308998 | 433 | -8 | A | America/Los_Angeles | airport | OurAirports | 0.027453 | 0.048108 |
| 3642 | 3830 | Chicago O'Hare International Airport | Chicago | United States | ORD | KORD | 41.978600 | -87.904800 | 672 | -6 | A | America/Chicago | airport | OurAirports | 0.060164 | 0.047940 |
| 191 | 193 | Lester B. Pearson International Airport | Toronto | Canada | YYZ | CYYZ | 43.677200 | -79.630600 | 569 | -5 | A | America/Toronto | airport | OurAirports | 0.042932 | 0.042651 |
# 例によって緯度・経度のデータ(pos)がある空港データに絞る
G_plot = G.subgraph(pos.keys())
# centralityの値を描画するノードのサイズに利用するために抽出。 1000は適当(描画の見栄えがよくなるように)
node_size_list = [node_data["betweenness_centrality"] * 1000 for node_name, node_data in G_plot.nodes(data=True)]
# 描画
fig = plt.figure(figsize=(15,7))
ax = fig.add_subplot(111)
nx.draw(G_plot, pos, ax, node_size=node_size_list, style=":")
plt.show()
次数中心性とは異なる結果が得られた。
例えばアラスカの真ん中あたりにある空港は、次数中心性の図では小さかったものの媒介中心性の図では大きい。この空港は規模こそ大きくないものの、航空網のつながりという観点では重要な空港である、といえるかもしれない。
(%%time をセルの一番上に書くと、実行時間を計測できる。)
%%time
clusters = greedy_modularity_communities(G)
CPU times: user 4.54 s, sys: 12.9 ms, total: 4.56 s Wall time: 4.55 s
# 中身の確認
for cluster_id, node_set in enumerate(clusters):
print(cluster_id, list(node_set)[:5])
0 ['KHH', 'PYY', 'LBS', 'SUG', 'LST'] 1 ['YVR', 'OCC', 'MSL', 'SHV', 'CYB'] 2 ['MHQ', 'WAW', 'TGP', 'TFS', 'ISU'] 3 ['HCR', 'NUP', 'VDZ', 'AKP', 'MNT'] 4 ['CBT', 'TNR', 'BKY', 'LUD', 'KAA'] 5 ['YPY', 'YYY', 'YBC', 'YCB', 'YEV'] 6 ['TQI', 'TQA', 'JUK', 'JUU', 'QOQ'] 7 ['SCZ', 'IRA', 'FRE', 'DLY', 'TAH'] 8 ['YNO', 'YTS', 'KEW', 'YAX', 'YHD'] 9 ['UAP', 'IPC', 'XMH', 'RKA', 'MAU'] 10 ['BYC', 'CIJ', 'UYU', 'JUL', 'RBQ'] 11 ['BUA', 'VAI', 'LSA', 'KVG', 'RAB'] 12 ['CKH', 'ODO', 'BQS', 'MJZ', 'YKS'] 13 ['CHX', 'BOC', 'SIC', 'PYC', 'OGM'] 14 ['SXK', 'OKL', 'SOQ', 'BIK', 'MKQ'] 15 ['NVK', 'MQN', 'SKN', 'BOO', 'SSJ'] 16 ['SVK', 'SQS', 'ORZ', 'SPR', 'TZA'] 17 ['HVG', 'SOJ', 'BVG', 'LKL', 'HFT'] 18 ['BMY', 'LIF', 'UVE', 'ILP', 'TGJ'] 19 ['NNM', 'CSH', 'VKT', 'UTS', 'KVX'] 20 ['ELG', 'TMR', 'OGX', 'DJG', 'INZ'] 21 ['FBS', 'LPS', 'YWH', 'RCE', 'CXH'] 22 ['EOI', 'NDY', 'NRL', 'WRY', 'SOY'] 23 ['ITB', 'MTE', 'STM', 'ORX', 'MEU'] 24 ['HOR', 'CVU', 'GRW', 'FLW', 'TER'] 25 ['INU', 'KSA', 'MAJ', 'TRW', 'PNI'] 26 ['TCD', 'LPD', 'ACR', 'LET', 'SVI'] 27 ['YDP', 'YNP', 'YHO', 'YRG', 'YMN'] 28 ['ERN', 'SJL', 'OLC', 'TFF', 'IRZ'] 29 ['WNR', 'ULP', 'BVI', 'BQL', 'BEU'] 30 ['KQA', 'IKO', 'AKB', 'DUT'] 31 ['ESD', 'BFI', 'CLM', 'FRD'] 32 ['BUC', 'ONG', 'DMD', 'NTN'] 33 ['AZI', 'ZDY', 'FJR', 'XSB'] 34 ['STZ', 'SXO', 'MQH', 'GRP'] 35 ['OND', 'ERS', 'MPA', 'NDU'] 36 ['TVY', 'MGZ', 'KAW'] 37 ['HFA', 'ETH', 'SDV'] 38 ['CMA', 'SGO', 'XTG'] 39 ['CMP', 'CDJ', 'RDC'] 40 ['PTJ', 'MEB', 'FLS'] 41 ['MNG', 'MGT', 'ELC'] 42 ['BTC', 'GIU', 'TRR'] 43 ['AJR', 'LYC'] 44 ['VHM', 'HMV'] 45 ['SSB', 'SPB'] 46 ['GCW', 'BLD'] 47 ['ABS', 'ASW'] 48 ['CKX', 'TKJ']
# clusters リストの内容を、 {(ノード名): (クラスタid)} の形式になるよう変換
node_cluster_dict = {}
for cluster_id, node_list in enumerate(clusters):
for node_name in node_list:
node_cluster_dict[node_name] = cluster_id
# 各ノードに cluster情報を付与
nx.set_node_attributes(G, node_cluster_dict, "cluster")
# 例によって緯度・経度のデータ(pos)がある空港データに絞る
G_plot = G_lcc.subgraph(pos.keys())
# 各クラスタに別々の色を付けたリスト(node_color_list)を作成 : わかりやすく色分けするため
color_list = list(mcolors.CSS4_COLORS.values()) # matplotlibの色一覧から色を抽出して、各クラスタに付与する色を取得
node_color_list = [color_list[node_data["cluster"]] for node_name, node_data in G_plot.nodes(data=True)]
# 各クラスタの中心座標を求める(ラベル付けのため)
cluster_pos = []
for cluster_nodes in clusters:
cluster_pos.append(airports[airports["IATA"].isin(cluster_nodes)][["Longitude", "Latitude"]].mean().values)
# 描画
fig = plt.figure(figsize=(15,7))
ax = fig.add_subplot(111)
## ネットワークの描画
nx.draw(G_plot, pos, ax, node_size=30, node_color=node_color_list, style=":")
## クラスタ番号の描画
for cluster_id, cluster_xy in enumerate(cluster_pos):
ax.annotate(
cluster_id,
xy=cluster_xy,
fontsize=10,
color='k',
bbox={'facecolor': color_list[cluster_id], 'edgecolor':'k', 'alpha':0.8}
)
plt.show()
https://sites.google.com/site/kztakemoto/r-seminar-on-igraph---supplementary-information#cartography より引用。
大規模なネットワークにおいてノードの役割や性質が分類できると便利ですよね。ひとつの考え方としてFunctional Cartography(Guimera & Amaral, 2005)という分類法があります。これはコミュニティ検出アルゴリズムなどをもちいてネットワークをモジュール分割し、モジュール内・外に対してノードがどのような役割を果たしているのかを特徴付けます。具体的には以下のように分類されます。 - R1: Ultra-peripheral nodes. 「超」周辺ノード - R2: Peripheral nodes. 周辺ノード - R3: Non-hub connectors. ハブではないが、複数のモジュールをつなぐコネクター - R4: Non-hub kinless nodes. ハブではないが、多くのモジュールにつながるノード - R5: Provincial hubs. モジュール内のハブ - R6: Connector hubs. 複数のモジュールをつなぐハブ - R7: Kinless hubs. 多くのモジュールにつながるハブ 例えば、ハブの存在するネットワークではハブの攻撃に対して脆弱(Albert et al., 2000)ですが、モジュール内のハブを攻撃するのと複数のモジュールをつなぐハブを攻撃するのとでは挙動が違います(Han et al., 2004)。このように、この分類で重要そうなノードを絞り込むことができます。
Functional Cartographyでは、コミュニティ検出によりモジュールに分けられた各ノードに対して、次の2つの値を求める。
ただし、
ただし、
これら2つの値を用いて、全ノードを以下のようにポジショニングし、分類する。
Functional Cartographyが提案された論文 : Functional cartography of complex metabolic networks
# 平均、標準偏差を求めるのでインポート
from statistics import mean, stdev
# 各ノードの Within-module degree を求める
within_module_degree_dict = {}
for cluster_nodes in clusters:
G_cluster_sub = G.subgraph(cluster_nodes)
cluster_degrees = [degree for node_name, degree in G_cluster_sub.degree()]
degree_mean = mean(cluster_degrees)
degree_stdev = stdev(cluster_degrees)
# 全ノードの次数が等しい場合には、標準偏差が0になり割り算でエラーが出るので、割らずに0を入れる
if degree_stdev == 0:
for node in cluster_nodes:
within_module_degree_dict[node] = 0
continue
for node in cluster_nodes:
within_module_degree_dict[node] = (G_cluster_sub.degree(node) - degree_mean) / degree_stdev
# Within-module degree 情報をネットワークに付与
nx.set_node_attributes(G, within_module_degree_dict, "within_module_degree")
# 各ノードの Participation coefficient を求める
participation_coefficient_dict = {}
for node in G.nodes:
node_p = 1
node_degree = G.degree(node)
node_neighbors = set(G.neighbors(node))
for cluster_nodes in clusters:
node_p -= (len(node_neighbors & cluster_nodes) / node_degree) ** 2
participation_coefficient_dict[node] = node_p
# Participation coefficient 情報をネットワークに付与
nx.set_node_attributes(G, participation_coefficient_dict, "participation_coefficient")
# 各ノードの役割(R1 ~ R7)を求める
node_part_dict = {}
for node_name, node_data in G.nodes(data=True):
if node_data["within_module_degree"] < 2.5:
if node_data["participation_coefficient"] < 0.05:
node_part_dict[node_name] = "R1"
elif node_data["participation_coefficient"] < 0.625:
node_part_dict[node_name] = "R2"
elif node_data["participation_coefficient"] < 0.8:
node_part_dict[node_name] = "R3"
else:
node_part_dict[node_name] = "R4"
else:
if node_data["participation_coefficient"] < 0.3:
node_part_dict[node_name] = "R5"
elif node_data["participation_coefficient"] < 0.75:
node_part_dict[node_name] = "R6"
else:
node_part_dict[node_name] = "R7"
# 各ノードの役割をネットワークに付与
nx.set_node_attributes(G, node_part_dict, "node_part")
# 各役割のノードがどのくらいあるか
Counter(node_part_dict.values())
Counter({'R1': 2831, 'R2': 457, 'R5': 85, 'R3': 16, 'R6': 36})
わかりやすくするため描画する
# ハブノードは色や大きさを変えて見やすくする
part_color_dict = {
"R1": "gray",
"R2": "blue",
"R3": "green",
"R5": "orange",
"R6": "red"
}
part_size_dict = {
"R1": 2,
"R2": 5,
"R3": 20,
"R5": 50,
"R6": 100
}
node_color_list = [part_color_dict[node_data["node_part"]] for node_name, node_data in G_plot.nodes(data=True)]
node_size_list = [part_size_dict[node_data["node_part"]] for node_name, node_data in G_plot.nodes(data=True)]
# 描画
fig = plt.figure(figsize=(15,7))
ax = fig.add_subplot(111)
## ネットワークの描画
nx.draw(G_plot, pos, ax, node_size=node_size_list, node_color=node_color_list, style=":")
plt.show()
Infomap法のpythonパッケージが公開されている( https://pypi.org/project/infomap/ )が、Windowsではパッケージをインストールできない。公式サイトにも
A pre-compiled version is available for macOS users.
と載っている。
そこで、今回はGoogle Colaboratoryを用いた方法を紹介する。 Google ColaboratoryではinfomapのPythonパッケージを使用することができる。
Google Colaboratoryの紹介 : https://gammasoft.jp/blog/google-colaboratory-for-learning/
routes.net) に出力する。# Infomapパッケージが読み込める形のファイルを出力。
nx.write_pajek(G, "routes.net")
続きはGoogle Colaboratoryでおこなう
https://colab.research.google.com/drive/1O1bPux0vjLyp7wvuDM4p86zhiLTnntCi?usp=sharing にコードがある
Google ColaboratoryでInfomapによるコミュニティ検出を行い、実行結果の modules.pickle ファイルが生成できたら、ダウンロードしてこのnotebookで読み込めるようにする。
# コミュニティ検出結果の読み込み
with open("modules.pickle", "rb") as f:
modules = pickle.load(f)
# 各ノードに付与するために{(ノード名): (クラスタid)} の形式になるよう変換
node_cluster_dict_infomap = {}
for node_name, module in zip(G.nodes, modules):
node_cluster_dict_infomap[node_name] = module
# 各ノードに cluster情報を付与
nx.set_node_attributes(G, node_cluster_dict_infomap, "cluster_infomap")
# 例によって緯度・経度のデータ(pos)がある空港データに絞る
G_plot = G.subgraph(pos.keys())
# 各クラスタに別々の色を付けたリスト(node_color_list)を作成 : わかりやすく色分けするため
color_list = list(mcolors.CSS4_COLORS.values())
node_color_list = [color_list[node_data["cluster_infomap"]] for node_name, node_data in G_plot.nodes(data=True)]
# 各コミュニティの中心座標を求める(ラベル付けのため)
cluster_pos = []
# コミュニティが13個あることを確認
Counter(modules)
Counter({1: 1993,
3: 218,
2: 914,
4: 183,
6: 37,
11: 2,
5: 52,
10: 4,
12: 2,
8: 4,
13: 2,
9: 4,
7: 10})
# 13個の各コミュニティの中心座標を集計
for cluster_id in range(1, 14):
node_cluster_filtered = list(filter(lambda x: x[1] == cluster_id, node_cluster_dict_infomap.items()))
cluster_nodes = set(map(lambda x: x[0], node_cluster_filtered))
cluster_pos.append(airports[airports["IATA"].isin(cluster_nodes)][["Longitude", "Latitude"]].mean().values)
# 描画
fig = plt.figure(figsize=(15,7))
ax = fig.add_subplot(111)
## ネットワークの描画
nx.draw(G_plot, pos, ax, node_size=30, node_color=node_color_list, style=":")
## クラスタ番号の描画
for i, cluster_xy in enumerate(cluster_pos):
cluster_id = i + 1 # クラスタ番号が1始まりになっているため、合わせる
ax.annotate(
cluster_id,
xy=cluster_xy,
fontsize=10,
color='k',
bbox={'facecolor': color_list[cluster_id], 'edgecolor':'k', 'alpha':0.8}
)
plt.show()
Modularity最大化とは違う分け方のコミュニティが得られた。